#张晗——big work
#安装包
BiocManager::install("DESeq2")
install.packages("ggrepel")
# 加载包
library(ggrepel)
library(DESeq2)   
library(ggplot2)
library("openxlsx")
#设置新路径
setwd("C:/Users/86183/Desktop/data for big work") 
#读取文件并将第一列设置为行名
data<-read.xlsx("GSE224239_genes_fpkm_expression.xlsx",rowNames = TRUE)

#直方图
#删除fpkm和其他列，保留count原始数据值
data1<-data[,-(1:14)]
data2<-unlist(data1)
#去除值为零的数据
data3<-data2[data2!=0]  
#取对数
data4<-log10(data3)
#画直方图
pdf("Histogram.pdf")
hist(data4,label=T,col = c("#0000FF","#FF1493"),
            main = "Histogram of GSE224239",xlab="data")
dev.off()    


countData <- data1
#获得countData数据
condition<-substring(colnames(countData),7,8)

#样本层次聚类
hdata=t(countData)
#计算行之间的距离
distance=dist(hdata, method = "manhattan") #曼哈顿距离计算样本间距离
tree<-hclust(distance,method="average")
pdf("hclust.pdf")
plot(tree,hang=-1,cex=.8)#下方对齐，点的大小为cex.8
dev.off()

#建立colData数据
condition<-factor(condition,labels = c("实验组","对照组"))
colData <- data.frame(row.names=colnames(countData), condition)
#调用包，差异分析
dds <- DESeqDataSetFromMatrix(countData = countData, 
                              colData = colData, 
                              design = ~ condition)
## 过滤低表达基因,至少有三个样本表达量高于10
dds <- dds[rowSums(counts(dds)>10)>=3,]
dds2<-DESeq(dds)
#用result()函数获取结果
res<-results(dds2)
#查看结果
summary(res)
# res格式转化：用data.frame转化为表格形式
res1 <- data.frame(res)
# 表达量显著上调的基因
res1_up<- res1[which(res1$log2FoldChange >= 1 & res1$pvalue < 0.05),]    
# 表达量显著下调的基因
res1_down<- res1[which(res1$log2FoldChange <= -1 & res1$pvalue < 0.05),]  
res1_total <- rbind(res1_up,res1_down)
#将上下调基因信息写入文件
write.csv(res1_up,file="DESeq2_results_up.csv")
write.csv(res1_down,file="DESeq2_results_down.csv")


genes<- res1
genes<-subset(genes,!is.na(genes$log2FoldChange)& !is.na(genes$padj))
# 根据上调、下调、不变为基因添加颜色信息
genes$group <- ifelse(genes$padj<0.05 & genes$log2FoldChange > 1,'up',
                    ifelse(genes$padj<0.05 & genes$log2FoldChange < -1,'down','normal'))
#将差异表达结果写入文件
write.csv(genes,file= "DESeq2_diffExpression.csv")

#火山图可视化
genes$label=rownames(gene)
genes$label=genename
pdf("volcano_plot.pdf")
p <- ggplot(
  # 指定数据、映射、颜色
  genes, aes(log2FoldChange, -log10(padj), col = group)) +  
  geom_point(size=2) +
  #重置颜色
  scale_color_manual(values = c(up="#FF1493",down = "#0000FF",normal = "darkgrey")) +
  # 标题、坐标轴名称
  labs(title="Volcano_plot of GSE224239",
       xlab="log2 (fold change)",ylab="-log10 (q-value)") +
  #画水平线
  geom_hline(yintercept = -log10(0.05), lty=4,lwd=1) +
  #画竖直线
  geom_vline(xintercept = c(-1, 1), lty=4,lwd=1) +
  #添加基因标签
  geom_text_repel(aes(label=label),size=2)+
  theme_bw()+
  #设置图标标题居中，设置图例靠右且不显示图例标题
  theme(plot.title=element_text(hjust=0.5), legend.position="right", 
        legend.title = element_blank())

p
dev.off()


#热图
library(pheatmap)
# 依次按照padj值log2FoldChange值进行排序
res1_total <- res1_total[order(res1_total$padj, res1_total$log2FoldChange, 
                               decreasing = c(FALSE, TRUE)), ]
#选择前四十个绘制热图
res_total=res1_total[1:40,]
df<-countData[intersect(rownames(countData),rownames(res_total)),]
df<-log2(df+1)
pdf("pheatmap.pdf")
pheatmap(df,
         #显示行名列名
         show_rownames = T,
         show_colnames =T,
         cluster_cols = T,#对列进行集群分析
         cluster_rows=F,#不对行进行集群分析
         height=13,  #高度为13
         main="pheatmap of GSE224239 ",#主标题
         frontsize = 9,#字体大小为9
         angle_col=90, #倾斜角度为90
         color =colorRampPalette(c("navy", "white", "firebrick3"))(100),
         clustering_method = 'single'#聚类方法
) 
dev.off()


#富集分析
if (!require("BiocManager",quietly=TRUE))
  install.packages("BiocManager")
if(!require("clusterProfiler",quietly = TRUE))
  BiocManager::install("clusterProfiler")
#安装大鼠基因库
BiocManager::install('org.Mm.eg.db')
library(org.Hs.eg.db)
library(clusterProfiler)
library(org.Mm.eg.db)
genes$ID=rownames(genes)
#筛选差异基因
ID=genes$ID[genes$padj<0.05]
ego=enrichGO(ID,"org.Mm.eg.db",
             qvalueCutoff = 1,
             pvalueCutoff = 1,
             pAdjustMethod = "none",
             keyType =  "ENSEMBL",
             ont="ALL",
             universe = genes$ID)
#将结果写入文件
write.table(ego,file="enrichGO.txt")
write.csv(ego,file="enrichGO.csv")

#结果可视化
#柱状图
pdf(file="eGO_barplot.pdf",width = 8,height = 10) 
barplot(ego, x = "GeneRatio", color = "p.adjust", 
        #默认参数（x和color可以根据eG里面的内容更改）
        showCategory =10, #只显示前10
        split="ONTOLOGY") + #以ONTOLOGY类型分开
  facet_grid(ONTOLOGY~., scale='free') #以ONTOLOGY类型分开绘图
dev.off()
#气泡图
pdf(file="eGO_dotplot.pdf",width = 8,height = 10) 
dotplot(ego,x = "GeneRatio", color = "p.adjust", size = "Count", #默认参数
        showCategory =8,#只显示前5
        split="ONTOLOGY") + #以ONTOLOGY类型分开
  facet_grid(ONTOLOGY~., scale='free') #以ONTOLOGY类型分屏绘图
dev.off()
#cnetplot
pdf(file="eGo_cnetplot.pdf",width = 8,height = 10) 
cnetplot(ego, showCategort=10,categorySize="pvalue",colorEdge = TRUE)
dev.off()
#emapplot
pdf(file="eGo_emapplot.pdf") 
compare_cluster_GO_emap <- enrichplot::pairwise_termsim(ego, semData = d)
emapplot(compare_cluster_GO_emap,showCategory = 20)
dev.off()
